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SUMMARY 

Coupled problems with various combinations of multiple physics, scales, and domains are found in 
numerous areas of science and engineering. A key challenge in the formulation and implementation 
of corresponding coupled numerical models is to facilitate the communication of information across 
physics, scale, and domain interfaces, as well as between the iterations of solvers used for response 
computations. In a probabilistic context, any information that is to be communicated between 
subproblems or iterations should be characterized by an appropriate probabilistic representation. 
Although the number of sources of uncertainty can be expected to be large in most coupled problems, 
our contention is that exchanged probabilistic information often resides in a considerably lower 
dimensional space than the sources themselves. In this work, we thus use a dimension-reduction 
technique for obtaining the representation of the exchanged information. The main subject of this 
work is the investigation of a measure-transformation technique that allows implementations to exploit 
this dimension reduction to achieve computational gains. The effectiveness of the proposed dimension- 
reduction and measure-transformation methodology is demonstrated through a multiphysics problem 
relevant to nuclear engineering. Copyright © 2000 John Wiley & Sons, Ltd. 
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1. Introduction 

The modeling and simulation of coupled systems governed by multiple physical processes that 
may exist simultaneously across multiple scales and domains are critical tools for addressing 
numerous challenges encountered in many areas of science and engineering. However, models 
are, by definition, only approximations of their target scenarios and are thus prone to 
modeling errors. Additionally, parametric uncertainties may exist owing to various limitations 
in manufacturing and experimental methods. Uncertainty quantification (UQ) thus constitutes 
a key requirement for achieving realistic predictive simulations. 
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Probability theory provides a rigorous mathematical framework for UQ, which permits 
a unified treatment of modeling errors and parametric uncertainties. The first step in a 
probabilistic UQ analysis typically involves using methods from mathematical statistics [1, 2] 
to characterize the uncertain features associated with a model as one or more random 
variables, random fields, random matrices, or random operators. The second step is to map 
this probabilistic representation of inputs through the system model into a probabilistic 
representation of responses. This can be achieved in several ways, which include Monte Carlo 
sampling techniques [3] and stochastic expansion methods. The latter typically involve the 
computation of a representation of the predictions as a polynomial chaos (PC) expansion. 
Several approaches are available to calculate the coefficients in this expansion, such as 
embedded projection [4, 5], nonintrusive projection [5], and collocation [6-10]. 

A key challenge in the formulation and implementation of a coupled model is to facilitate the 
communication of information across physics, scale, and domain interfaces, as well as between 
the iterations of solvers used for response computations. This information can comprise physical 
properties, energetic quantities, or solution patches, among other quantities. Although the 
number of sources of uncertainty can be expected to be large in most coupled problems, we 
believe that the exchanged information often resides in a considerably lower dimensional space 
than the sources themselves. Exchanged information can be expected to have a low effective 
stochastic dimension in multiphysics problems when this information consists of a solution 
field that has been smoothed by a forward operator and in multiscale problems when this 
information is obtained by summarizing fine-scale quantities into a coarse-scale representation. 

In a previous paper [11], we had proposed the use of a dimension-reduction technique to 
represent the exchanged information: we proposed to represent the exchanged information by 
an adaptation of the Karhunen-Loeve (KL) decomposition as this information passes from 
subproblem to subproblem and from iteration to iteration. The main objective of this paper is 
to complement this dimension-reduction technique by a measure-transformation technique 
that allows implementations to exploit the dimension reduction to achieve computational 
gains. Here, we specifically consider implementations that use stochastic expansion methods. 
The proposed measure-transformation technique allows such implementations to carry out key 
algorithmic operations, such as the construction of PC expansions, with respect to the reduced 
stochastic degrees of freedom of the exchanged uncertainty representations, thus leading to a 
solution in a reduced-dimensional space, which in turn reduces the computational cost. 

The organization of this paper is as follows. First, in Sec. 2, we outline the proposed 
dimension-reduction and measure-transformation methodology. Then, in Sees. 3-5, we present 
the dimension-reduction and measure-transformation techniques. In Sec. 6, we provide details 
on the implementation of these techniques. Finally, in Sees. 7 and 8, we demonstrate the 
proposed methodology through an illustration problem. 



2. Dimension- reduction and measure-transformation methodology 
2.1. Model problem 

This paper is devoted to the solution of a stochastic coupled model of the following form: 
f(u, x, £) = 0, y = h(u, £), / : E r x R s " x K m -> R r , h : R r x K m -> R r ° , 
g{y,v,C) = Q, x = k(v,C), 9 '■ R r ° x f x P 4 I s , k : K s x E" -> R s °. 
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To avoid certain technicalities involved in infinite-dimensional representations, we assume that 
these equations are discretized representations of a stochastic model that couples two physics, 
two scales, two domains, or a combination of these sub-problems. For instance, these equations 
may be obtained from the spatial discretization of a steady-state problem, or they may be 
the equations obtained at a single time step after the spatial and temporal discretization of 
an evolution problem. Further, we assume that the data of the first subproblem, which enter 
this subproblem as coefficients or loadings or both, depend on a finite number of uncertain 
real parameters denoted as £1, . . . ,£ m , and that the data of the second subproblem depend 
on a finite number of uncertain real parameters denoted as d, . . . , £„. Lastly, we model these 
sources of uncertainty as random variables and collect them into vectors, £ = (£i,...,£m) 
and C — (Cij • • • ; Cn), which are assumed to be defined on a probability triple (0, T, P) and 
considered to have values in R m and R™, respectively. 

The stochastic coupled model (1) is a general bidirectionally coupled model. The solution 
variables, u, of the first subproblem, /, depend on the solution variables, v, of the second 
subproblem, g, through coupling variables, x; likewise, v depends on u through y. 

Thus, to solve this stochastic coupled model, we require to find the random variables u and v 
defined on (0, T, P) with values in W and R s such that (1) is satisfied under the assumption 
that the stochastic coupled model is well-posed in that it admits a unique and stable solution. 

2.2. Partitioned iterative solution 

Because a coupled model usually characterizes its response only in an implicit manner, the 
numerical solution of a coupled model typically requires an iterative method. Here, we assume 
that iterative methods and associated solvers already exist for each subproblem, and we 
therefore consider a partitioned iterative method that reuses the aforementioned separate 
solvers as steps in a global iterative method built around them to obtain a solution to the 
coupled model. Let us assume that each of the aforementioned separate iterative methods is 
based on the reformulation of the associated subproblem as a fixed-point problem: 

u = a(u,x,£), y = h{u,i), a : R r x R S( > x R m -> R r , h : W x R m -> R r °, 
v = b(y,v,C), x = k(v,C), b : R r ° x R s x R" -> R s , k : R s x R" -> R s ° . 

It should be noted that these equations can be obtained by setting a(u, v, £) = u — f(u, v, £) 
and b(u, v, £) = v — g(u, v, £), but that alternative reformulations, such as those involving a 
direct solution of the subproblcms or of their linear approximations, are often better adapted. 
We then consider the solution of the stochastic coupled model by a Gauss-Seidel iterative 
method using suitable initial values u°, v°, and x° = k(v°, £) as follows: 



This is not the only partitioned iterative method available; however, for simplicity, we employ 
only this method in this work. It should be noted that although we implement the proposed 
methodology using the Gauss-Seidel iterative method, one can readily use the proposed 
methodology with other iterative methods such as Jacobi, relaxation, and Newton methods. 




aiu'-W- 1 ,*), 



y l = h(u£,Z), 

x e = k(v e ,C). 



(3) 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



STOCHASTIC MODELING OF COUPLED PROBLEMS 



3 



2.3. Dimension reduction 

We believe that exchanged information often resides in a considerably lower dimensional 
space than the input sources of uncertainty themselves. Therefore, in [11], we investigated 
the effectiveness of dimension-reduction techniques for the representation of the exchanged 
information. Rather than exchanging the coupling variables x l and y l and the solution 
variables u and v in their original form, we proposed to approximate these random variables 
by a truncated KL decomposition as they pass from subproblcm to subproblcm and from 
iteration to iteration. The use of this dimension-reduction technique leads to the solution of 
the stochastic coupled model by a Gauss-Seidel iterative method as follows: 

u* = a(u*-^,x*-^,t), y l = h{u\i), 

where q l ' d = [y i,d ; v e ~ 1 > d ] and r l - e = [M £_1 ' e ; & £_1,e ] are truncated KL decompositions of 
q l = [y 1 -.^- 1 ] and r l — respectively, which are written as follows: 

T 1 (5) 

These decompositions are described in detail in later sections. Truncation of KL decompositions 
most often results in approximation errors, and we thus use a hat superscript to distinguish 
the successive approximations determined by (4) from those determined by (3). 

It should be noted that (5) provides a combined reduced-dimensional representation of y E 
and v e ~ 1 in terms of a single set of reduced random variables rf — {r}{, . . . , r] d ) and a combined 
reduced-dimensional representation of i/^ 1 and x^~ x in terms of a single set of reduced random 
variables </ = (t[, . . . , i*). However, this is not the only construction of a reduced-dimensional 
representation that could be considered. One can readily use the proposed methodology 
with other dimension-reduction techniques such as those involving the construction of a 
separate reduced-dimensional representation of the coupling and solution variables, with each 
representation having its own reduced random variables. 

Finally, although our notations do not express a potential dependence of d and e on £, it 
should be noted that the reduced dimensions can be allowed to depend on the iteration. 

2.4- Measure transformation 

The successive approximations determined by the iterative method (3) that does not involve 
dimension reduction can be constructed as random variables of the following form: 

u<(0) = u<(*(0),C(0)), 

i.e., u and v e can be constructed as transformations of the input random variables £ = 
(£i, . . . , £ m ) and £ = (£i, . . . , ( n ). The random variables u l and v l thus exist in a solution space 
of stochastic dimension m+n. Consequently, implementations of (3) using stochastic expansion 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



4 



M. ARNST, R. GHANEM, E. PHIPPS AND J. RED-HORSE 



methods would typically involve the approximation of u l and v l by finite-dimensional 
representations as follows: 



< e R\ 

(7) 

vt g R s , 



\ a \=0 
|c|=0 

where {ip a ,a e i s a Hilbertian basis for the Hilbert space of Pr^ ^)-square-integrable 

functions from R m+ ™ into R, in which P(£X) denotes the joint probability distribution of the 
input random variables £ and £, and |cc| = a± + . . . + a m+n . Then, the task of the solution 
algorithm would be to compute the coordinates u l a and v l a in these expansions. Depending 
on the specific nonintrusive projection, embedded projection, or collocation method that is 
chosen, this task would require the construction of basis functions, quadrature rules, moment 
tensors, or a combination of these with respect to the joint probability distribution P(£X)- 

In contrast, owing to the dimension reduction, the successive approximations determined by 
the iterative method (4) can be constructed as random variables of the following form: 

«'(0) = u'(««) )l '((l)), 

b'(0) = 6V(0),C(0)); 
i.e., ii l can be constructed as a transformation of the input random variables £ = (£i, . . . , £ m ) 
and the reduced random variables </ = and v l can be constructed as a 

transformation of rf — (r][ , . . . , 77^) and C — (Cii • • ■ > Cn)- Thus, the random variable ii l exists 
in a space of stochastic dimension m + e and the random variable v l exists in a space of 
stochastic dimension d+n. Consequently, implementations can exploit the dimension reduction 
to approximate u l and v l by finite-dimensional representations as follows: 



(9) 



|/3|=0 

ItI=o 

where {I^,/3 e N m+e } and e N d+ ™} are Hilbertian bases for the Hilbert spaces 

of P(^ jt i)-square-integrable functions from M m+e and of ^-square integrable functions 
from R d+ ™, respectively, into R. Here, P(£, t t) and P( v e x) are the joint probability distributions 
of £ and l 1 and of rf and £, respectively. Now, depending on the specific solution method that 
is chosen, the construction of these expansions requires the construction of basis functions, 
quadrature rules, moment tensors, or a combination of these with respect to the joint 
probability distributions P(^^) and P^ q. This approach involves a measure transformation 
because from the dimension reduction given by (5), it follows that the probability distributions 
of the reduced random variables i l and rf are transformations of the probability distribution 
of the input random variables £ and £, as described in detail in a later section. 

Throughout this work, we have employed Hilbertian bases that are constituted of 
polynomials and we thus refer to them as polynomial chaos (PC) bases. 
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Although our notations do not express a potential dependence of p and q on the subproblcm 
or £, it should be noted that the subsets of basis functions used to construct the finite- 
dimensional representations can be allowed to depend on the subproblem and the iteration. 

2.5. Effectiveness of the dimension-reduction and measure-transformation methodology 

The key feature of the proposed methodology is that it enables a solution of the subproblems 
in a reduced-dimensional space when the exchanged information has a low effective stochastic 
dimension. Specifically, a solution in a reduced-dimensional space is enabled when the reduced 
dimensions can be selected such that d < m and e < n while sufficient accuracy is maintained; 
refer to (6) and (8). This benefit is of particular significance for implementations of stochastic 
coupled models that use stochastic expansion methods. These methods suffer from a curse 
of dimensionality in that their computational cost increases quickly with an increase in the 
stochastic dimension. The proposed methodology addresses the curse of dimensionality by 
mitigating the increase in stochastic dimension when information is exchanged. 

Finally, it should be noted that the proposed methodology can readily be adapted to meet 
various requirements of specific applications. One could, for instance, use the KL decomposition 
to represent only those exchanged random variables that are of low effective stochastic 
dimension and use another representation for the remaining variables; further, one could 
implement a measure transformation to solve only those subproblems whose computational 
cost would thus be lowered and implement another approach for the remaining subproblems. 



3. Karhunen-Loeve decomposition 

Here, we concisely recall the KL decomposition used in [f f] to construct a reduced-dimensional 
representation of a random variable q that is defined on a probability triple (@,T,P), takes 
values in a finite-dimensional Euclidean space M. w , and is of the second order: 

/ ||<z|| 2 dP<+oo, (fO) 

where ||-|| denotes the Euclidean norm. Refer to [f f] and the references therein for more details 
on this construction. 

3.1. Second-order descriptors 

The mean vector q and covariance matrix C q of the second-order random variable q are defined 
as the w-dimensional vector and square matrix, respectively, such that 

q= f qdP, (11) 
C q = f {q-q){q-q) T dP. (12) 

Je 

3.2. Reduced- dimensional representation 

Let W be a w-dimensional, square, symmetric, and positive definite matrix; we refer to this 
matrix as the weighting matrix and comment on its usefulness later. Then, because C q is 
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symmetric and positive scmidcfinite, the solution of the generalized eigenproblcm 

W T C q W<l>> = \jW4j (13) 

provides a set of w eigenvalues Ai > A2 > ... > \ w > in addition to w 
eigenvectors 1 , . . . , <p w , which constitute a W- weighted orthonormal basis of W" such that 

(^) T W=<5y, (14) 

where 5ij is the Kronecker delta, which is equal to 1 if i = j and otherwise. The KL-type 
decomposition of q is then given by 



q = q + ^2^Vj<l> 3 , (15) 

where the rjj are random variables defined on (9, T, P), with values in R, such that 

r h = -^=(q-qfW&. (16) 
V A i 



These reduced random variables r]j are zero-mean and uncorrelated: 

/ VjdP = 0, (17) 
Je 

[ r ll r lj dP = S lJ . (18) 
The truncation of (15) after d terms provides a reduced-dimensional representation as follows: 

d 

Although the reduced random variables of the KL decomposition are uncorrelated, it should be 
noted that they are in general statistically dependent. Further, although the joint probability 
distribution of the reduced random variables is usually complicated, it is determined by the 
KL decomposition. Finally, if the random variable to be reduced is Gaussian, then the reduced 
random variables are also Gaussian; however, the joint probability distribution of the reduced 
random variables is most often not a "labeled" probability distribution. 



3.3. Concluding remarks 

In this section, we presented a version of the KL decomposition that features a weighting 
matrix. We had demonstrated in [11] that this weighting matrix is particularly useful for the 
construction of a reduced-dimensional representation of a random variable that solves a space- 
time discretized stochastic model: we showed that by appropriately choosing the weighting 
matrix as the Gram matrix of the discretization basis, a reduced-dimensional representation 
is obtained which is consistent with the function-analytic structure that the stochastic model 
exhibited prior to its discretization. Nevertheless, it should be noted that the standard KL 
decomposition is recovered by simply setting the weighting matrix equal to the identity matrix. 
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4. Polynomial chaos with respect to arbitrary probability distributions 

Here, we consider the construction of a PC basis with respect to an arbitrary, although fully 
specified, probability distribution P x defined on a finite-dimensional Euclidean space W . 

4-1. Objective 

In the context of the proposed methodology, wc consider P x as the probability distribution 
of a random variable x that collects the sources of uncertainty that enter a subproblem of 
a stochastic coupled model. With reference to (9), at a specific iteration, x could collect the 
components of £ and l with z = m + e or the components of 77 and £ with z = d + n, in which 
case the sought PC basis would be needed to build a PC expansion of u or v, respectively. 

4.2. Methodology 

If the probability distribution P x does not exhibit statistical dependence and its marginal 
probability distributions are "labeled" univariate probability distributions, then a PC basis 
can readily be constructed by tensorization of the classical univariate orthonormal polynomials 
that are associated with these univariate probability distributions in the Askey table [4, 5,12- 
14] . Well-known examples include the use of PC bases of Hermite and Legendre polynomials 
with respect to Gaussian and uniform probability distributions, respectively. 

However, here, we are interested in a general case wherein the probability distribution P x 
may exhibit statistical dependence and may not be "labeled." This case is considered in 
the context of the proposed methodology because the reduced random variables of a KL 
decomposition are in general statistically dependent and not "labeled." 

If P x exhibits statistical dependence or is not "labeled," the Hilbert space of P x -square- 
intcgrable functions from M z into K does not necessarily admit a PC basis, as discussed in 
the following sections. Further, polynomials that are orthonormal with respect to an arbitrary 
probability distribution usually cannot be read from tables in the literature, thus requiring 
a computational construction. Therefore, the approach adopted in this section is as follows. 
First, we describe a computational construction of a set of orthonormal polynomials. Then, 
we discuss conditions under which this set is also a PC basis. 

4-3. Set constituted of P x - orthonormal polynomials 

We use standard notations involving multi- indices to work with multivariate polynomials. 
We refer to elements (3 — . . . ,(3 Z ) of N z as multi-indices. A (multivariate) monomial x^ 
associated with a multi-index (3 is then a function from M. z into M defined by 

X P = X f t x...xxf*. (20) 

Further, we refer to the number |/3| = /3i + . . . + /3 Z as the modulus of the multi-index (3 and 
also as the total degree of the corresponding monomial x P ■ A (multivariate) polynomial is then 
a function from M z into R that maps any x to a finite sum cpx with real coefficients cp. 
Let V% be the space of all polynomials in z variables with a total degree of at most q: 

VI = U : X = 08i,---,&)->t(x)= ]T c pX P \ (21) 

L l/3|=0 ' 
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the dimension of V% is commonly expressed as 

dto W )-(* + «)-<f + f>!. (22) 

In order to construct a set of P x -orthonormal polynomials in V% for a given total degree q, we 
propose a procedure wherein we first arrange the monomials {x 13 ', < |/3| < q} spanning V% 
in a sequence and then orthonormalize this sequence by the Gram-Schmidt method as follows: 

ffl(x) = /T%f 7 (|3), with T Z ix)dPx , (23) 

Tp(x) = fp(x)/ I f p {x)fp{ X )dP x . (24) 

Clearly, the Gram-Schmidt method cannot always be executed without difficulty: it will 
fail when the integral in the numerator of (23) does not exist or when the integral in the 
denominator of (24) vanishes. The former difficulty may occur when a few or all the moments 
of P x are unbounded. The latter difficulty corresponds to the case wherein Tp has a vanishing 
weight with respect to P x and may occur either when P x is a discrete probability distribution 
or when P x is degenerate in that it is concentrated on a hypersurface in R z . 

One can readily impose suitable conditions on P x to ensure proper execution of the Gram- 
Schmidt method. Indeed, if P x has finite moments up to a total degree of 2q, i.e., 

/ Wi x...x X f*|dP x <+oo, 0<|/3|<2g (25) 

and if P x is nondegenerate, i.e., 

/ 7r(x) 2 dP x > 0, ViePf, TT^O, (26) 

then the aforementioned difficulties cannot occur; thus, it is then ensured that the Gram- 
Schmidt method provides a unique set of polynomials Tp that are mutually orthogonal in 
that J RZ Tp(x)T-y(x)dP x = if f3 ^ 7 and normalized so that f RZ \Tp(x)\ 2 dP x = 1. 

A frequently encountered example of a case wherein P x has finite moments and is non- 
degenerate is the one in which P x admits a probability density function that has a closed and 
bounded support with a nonempty interior. 

It should be noted that cases wherein P x is a discrete or degenerate probability distribution 
can easily be handled by simply discarding polynomials with a vanishing P x -weighted norm. 

4-4- PC basis constituted of P x -orthonormal polynomials 

The probability distribution P x determines a space L 2 p (M z , R) of functions from M z into R 
that are square- integrable with respect to P x : 

f \f(x)\ 2 dP x <+^ V/6Lj, (R',R). (27) 

The function space L 2 p ^ (R z , R) is a Hilbert space for the inner product 

(/,<?>= / f(x)9(x)dP x , Vf,geL Px (R z ,R). (28) 
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If the conditions (25) and (26) are fulfilled for a specific total degree q, the Gram-Schmidt 
method given by (23) and (24) provides a unique set {Tp,0 < |/3| < q} of P x -orthonormal 
polynomials spanning V%. This set of polynomials allows one to associate to each function / 
in L 2 p ^(R z ,R) a corresponding set {fp,0 < |/3| < q} of coordinates fp in R defined by 



U = (/.^) = / f(x)Tp(x)dP x , < |/3| < q. (29) 
This set of coordinates in turn determines an approximation f q of / in V% as follows: 

/«/*= E Z^. (30) 
101=0 

Clearly, if the conditions (25) and (26) are fulfilled for any total degree q, then the above 
construction provides a q- indexed family of the approximations f q of /. Then , the convergence 
of these approximations with respect to the total degree q is a critical issue. Ideally, one would 
like the approximation f q to converge to / as q is increased, i.e., 



lim / f(x) - E ffPM 



dP x = 0, V/e4jt 2 ,R). (31) 

1/31=0 

Now, a set {T/3, /3 € N z } of P x -orthonormal functions in Lp (R Z ,R) (which need not 
necessarily be polynomials) is called [15] a Hilbertian basis for Lp^(R z ,M) if it satisfies 

/ = E V/ei£ x (R*,R) (32) 

in the sense that the right-hand side should converge to the left-hand side in the L 2 -norm 
without depending on the ordering of the mult i- indices. Clearly, the fulfillment of (32) for 
any ordering implies the fulfillment of (31) for that specific ordering by the total degree. As 
mentioned previously, we refer to a Hilbertian basis constituted of polynomials as a PC basis. 

The family of Hermite polynomials is well known to provide a PC basis for a Hilbert space of 
functions that arc square-integrable with respect to a Gaussian probability distribution [13]. In 
general, the families of classical orthonormal polynomials are well known to provide PC bases 
for Hilbert spaces of functions that are square-integrable with respect to the corresponding 
"labeled" probability distributions in the Askey table [4, 5, 13, 14]. However, it is not obvious 
for a set of orthonormal polynomials, such as the one determined by the Gram-Schmidt method 
given by (23) and (24), to be a PC basis; this property generally depends on the probability 
distribution under consideration. 

If the support of the probability distribution P x is a nonempty, closed, and bounded subset 
of the z-dimensional Euclidean space R z , then by the Stone- Weierstrass theorem [16], any set 
of P x -orthonormal polynomials that spans the space of all the z-variate polynomials, including 
the set of orthonormal polynomials determined by the Gram-Schmidt method given by (23) 
and (24), forms a PC basis for the Hilbert space of functions that are square-integrable with 
respect to this probability distribution. The illustration problem considered in the last two 
sections of this paper falls within this case. 

On the other hand, cases wherein P x is a probability distribution with an unbounded 
support are more complex. For brevity, here, we only give one example to illustrate that sets 
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of orthonormal polynomials denned with respect to a probability distribution with unbounded 
support need not form a PC basis: consider a Gaussian random variable £; it can then be 
shown [17, 18] that the Hilbert space of functions from M into R that are square-integrable 
with respect to the probability distribution of \ — £ 2fe+1 with k — 1,2... does not admit a 
PC basis. Refer to [19] and the references therein for further details on this issue. 

4-5. Bibliographical comments 

Orthogonal polynomials have been a subject of intensive research in the literature; for instance, 
refer to the recent references [20-22]. In the research area of stochastic modeling and analysis, 
the construction of PC bases with respect to arbitrary probability distributions has already 
been considered in [23-25]: reference [23] presents a multielement PC basis synthesized from 
univariate orthonormal polynomials generated by the Stieltjes method; the method given in [24] 
involves a multidimensional PC basis deduced by tensorization of univariate orthonormal 
polynomials generated by the Gram-Schmidt method; and reference [25] proposes the use 
of a singular value decomposition to obtain a multidimensional PC basis that is orthonormal 
with respect to a discrete approximation of the probability distribution. 



5. Quadrature rules with respect to arbitrary probability distributions 

In this section, we consider the construction of a family of quadrature rules for the evaluation 
of the integrals of continuous functions from a finite-dimensional Euclidean space R z into R 
with respect to an arbitrary, although fully specified, probability distribution P x on R z . 

5.1. Objective 

Again, we consider P x as the probability distribution of a random variable x that collects the 
sources of uncertainty that enter a subproblcm of a stochastic coupled model. With reference 
to (9), at a specific iteration, x could collect the components of £ and l with z = m + e or the 
components of rj and £ with z = d + n. The purpose of the family of quadrature rules depends 
on the specific stochastic expansion method that is selected for discretization. For instance, if a 
nonintrusive projection method were adopted, a family of quadrature rules would be required 
to evaluate the projection of the solution to the subproblem onto the employed PC basis. 

5.2. Methodology 

If the probability distribution P x does not exhibit statistical dependence, then an appropriate 
family of quadrature rules can often be constructed through a full or sparse-grid tensorization 
of the families of Gaussian quadrature rules associated with these univariate probability 
distributions. Gaussian quadrature rules can be read from tables in the literature for "labeled" 
univariate probability distributions and can be computed easily for arbitrary univariate 
probability distributions [22]. Well-known examples include the use of Gauss- Hermite 
and Gauss-Legendre quadrature rules with respect to Gaussian and uniform probability 
distributions, respectively. The families of quadrature rules thus obtained are well known to 
be efficient in that they provide fast convergence rates for smooth integrands. 

However, here, we are interested in a general case wherein the probability distribution P x 
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may exhibit statistical dependence. This case is considered here because the reduced random 
variables of a KL decomposition are generally statistically dependent. 

General multivariate probability distributions do not necessarily admit Gaussian quadrature 
rules, as discussed in the following sections. Further, quadrature rules for integration with 
respect to arbitrary probability distributions usually cannot be read from tables in the 
literature, thus requiring a computational construction. Thus, the approach adopted in this 
section is as follows. First, we propose a method for the computational construction of a family 
of quadrature rules for the evaluation of the integrals of continuous functions from R z into R 
with respect to the probability distribution P x . Then, we discuss the efficiency and convergence 
properties of this method in addition to its relationship to other available methods. 

5.3. Proposed computational construction of quadrature rules 

Let A > 1 be a fixed finite integer that determines the level of the quadrature rule to be 
constructed by requiring its polynomial degree of exactness to be 2 A — 1. Further, let P x be a 
probability distribution on M z that satisfies 

/ Ixf 1 x ... xx*'\dP x <+™, 0<|/3|<2A-1. (33) 

Jut* 

We then propose the following two-step strategy for the construction of a quadrature rule: 

• First, we construct a quadrature rule of the following form: 

V 

f(x)dP x nJ2f(Xk)tik, (34) 

which allows integrals relative to P x to be approximated accurately, but may have a very 
large number v of nodes X\i ■ ■ ■ > Xp an d associated weights w\, . . . , wp. 

• Then, we construct an embedded quadrature rule of the following form: 

f{x)dP x R>52f(Xk)«>k, (35) 
which has level A in that it integrates all polynomials up to total degree 2A — 1 exactly: 

/ 7r(x)dP* = £>(XfcW> VTrePf" 1 , (36) 

J *' fc=i 

but uses only a small subset Xi = Xki > ■ • ■ > Xv = Xk v of ^ <C £ nodes of the original rule. 

The nodes and weights of the original quadrature rule (34) could, for instance, be constructed 
by a Monte Carlo integration approach. Here, the key challenge is rather in the selection of 
the small subset of nodes to be retained by the embedded quadrature rule. Clearly, in order 
to solve this subset selection problem, we require to determine a set of weights zu\ , . . . , wp for 
the candidate nodes respectively, such that these weights are mostly zero while 

they still allow integrals of polynomials up to the prescribed total degree of 2 A — 1 to be 
evaluated exactly. We propose to construct such a sparse set of weights through the solution 
of the following L 1 -minimization problem: 

min ||ct|| l i , subject to Azu = b, (37) 



/ 

Jr 



/ 

Jr 
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where ||-Cc7|| L i = |n7i| + \vj2\ + . . . + |ra7j>_i| + \wo\ and the equality constraints Azu = b serve 
to ensure that polynomials up to the prescribed total degree of 2A — 1 are integrated exactly: 

: : : : (38) 

b=[f v *i(x)dP x ... J R ,^(x)dP x ] T . (39) 

Here, m , . . . , ir^ is a basis for Vl x ~ x with = dim^f A_1 ). Let Bl)J\f then be a partitioning of 
the index set {1, . . . , v) associated with a solution vj to (37) such that the fc-th component Wk 
of that solution is nonzero when k is included in B and vanishes when k is included in N. 
The nodes of the original quadrature rule labeled by the indices in B = {k\, . . . , k u } are then 
the nodes Xi — Xki > • • • > X v = Xk„ *° retained by the embedded quadrature rule with the 
associated weights w\ — zu^ , ■ . ■ , w v — zuk„ , while the nodes labeled by the indices in N are 
to be discarded. In Sec. 6, we have demonstrated how such a sparse solution can be obtained 
using either a simplex [26] or an interior-point optimization algorithm [26-28] . Both of these 
approaches yield a sparse solution that has at most dim('p2 A ~ 1 ) nonzero components. Thus, 
the proposed construction yields a quadrature rule that uses at most as many nodes as there 
are equality constraints imposed in (37) to ensure the prescribed polynomial exactness. 

5.4. Efficiency 

The proposed construction yields a quadrature rule that uses at most dim('Pf *) nodes to 
achieve exactness for all z-variate polynomials up to a total degree of 2 A — 1. An important 
question thus pertains to the optimality of the proposed construction with respect to the 
number of nodes used to achieve this degree of exactness. For quadrature rules used to evaluate 
one-dimensional integrals, it is well known that a A-point Gaussian quadrature rule has a 
degree of exactness of 2A — 1 and that there exists no A-point quadrature rule that is exact 
for all polynomials up to a degree of 2A. In contrast, for quadrature rules used to evaluate 
multidimensional integrals, the minimum number of nodes required for the exact integration 
of a given number of polynomials is yet unknown in current state-of-the-art mathematics. 

Nevertheless, if the probability distribution P x has a closed and bounded support, it is 
known that the minimum number of nodes required for a quadrature rule to have a degree 
of exactness of 2A — 1 is [29] greater than or equal to dim^^ 1 ) and less than or equal 
to dim(7 7 | A_1 ). This upper bound corresponds to the Tchakaloff theorem: 

Theorem 5.1. (Tchakaloff [SO]) Let P x be a probability distribution on W z that admits a 
probability density function that is supported by a closed and bounded subset JC o/K z . Then, 
for any integer 2A — 1 > 0, there exists fj, = dimifPi^' 1 ) nodes Xi, ■ ■ ■ > X^ in JC and positive 
weights w\, . . . , such that the resulting quadrature rule has a degree of exactness of 2 A — 1. 

Hence, our construction provides a quadrature rule that uses at most as many nodes as the 
Tchakaloff theorem indicates are at most required to achieve the prescribed degree of exactness. 
It should be noted that this theorem not only ensures the existence of a dim(V^ A_1 )-point 
quadrature rule with a degree of exactness of 2 A — 1, but also guarantees the existence of such 
a rule with only positive weights. In contrast, our construction does not necessarily yield a 
quadrature rule with only positive weights. Further, it should be noted that extensions of this 
theorem to cases involving probability distributions with unbounded support exist [31]. 
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5. 5. Convergence 

When P x has finite moments of any order and the assumption (25) is thus fulfilled for any 
level A, the proposed construction provides a A-indexed family of quadrature rules, each 
member of which is required to achieve a corresponding degree of exactness of 2 A — 1. An 
important question then pertains to the convergence of the approximation of an integral of a 
given function using a sequence of quadrature rules with an increasing degree of exactness. 

If the probability distribution P x has a closed and bounded support, an insightful, although 
not comprehensive, answer is contained in the following theorem: 

Theorem 5.2. (Convergence of polynomial-based quadrature rules [32]) Let P x be a 
probability distribution on M. z that admits a probability density function that is supported by a 
closed and bounded subset JC of W . Then, a quadrature rule with degree of exactness of 2A — 1 
satisfies for every continuous function f from JC C K z into R the following inequality: 

f f(x)dP x -iTf(x k )w k < (l + J2\w k \) mm_ i max\f{ X )-*(x)\- (40) 

The proof follows from the triangle inequality; in fact, the left-hand side satisfies the inequality 



I r c II/* v I u v 

1-h.s. < f(x)dP x - h( X )dP x \ + \ L(x)dP x -^2n( Xk )w k + /( Xfc ; 



(41) 



<max xeK |/( X )-^(x)| =0 <E^ =1 k fc |max xex; |/(x)-vr(x)| 



The result (40) shows that the quadrature error is bounded by a product of two factors; the 
first grows with the sum of the absolute values of the weights, and the second is a bound 
on the approximation error introduced by the best approximation of the integrand by a 
polynomial with a total degree of at most 2 A — 1. This result justifies our construction (37) 
of the embedded quadrature rule. Our choice of the objective function as the sum of the 
absolute values of the weights results in the minimization of the first factor mentioned above, 
contributing to the upper bound on the quadrature error. Our requirement for the embedded 
quadrature rule to have a polynomial degree of exactness of 2A — 1 precisely results in the 
second factor mentioned above. This requirement allows the embedded quadrature rule to 
achieve a fast convergence rate for smooth integrands as the polynomial degree of exactness 
increases because the approximation error introduced by the best approximation of a smooth 
function by a polynomial of a specified degree on a closed and bounded set is well known to 
decrease at a fast rate with an increase in this degree. 



5.6. Bibliographical comments 

Three classes of approaches for multivariate integration have received most attention in the 
literature, namely, probabilistic and number-theoretic methods, polynomial-based methods, 
and adaptive techniques. Probabilistic and number-theoretic methods include Monte Carlo 
and quasi- Monte Carlo integration; for instance, refer to [3, 33]. These methods are well 
known to exhibit a rather slow rate of convergence as a function of the number of integrand 
evaluations, but Monte Carlo methods have the advantage that their rate of convergence is 
independent of the stochastic dimension. Polynomial-based integration rules are designed to 
be exact for a prescribed collection of polynomials; for instance, refer to [20, 29, 32, 34]. 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



14 



M. ARNST, R. GHANEM, E. PHIPPS AND J. RED-HORSE 



Polynomial-based methods have the advantage that their rate of convergence increases rapidly 
with the smoothness of the integrand, but these methods suffer from a curse of dimensionality 
in that their rate of convergence decreases with an increase in the dimension. Adaptive methods 
involve the choice of the nodes and weights in a manner that is dependent on the integrand 
to achieve fast convergence rates, while still limiting the increase in computational cost as the 
dimension increases; for instance, refer to [35]. 

Among the polynomial-based methods, four classes of approaches have been investigated 
extensively in the literature. The first class includes methods that rely on the direct numerical 
solution of the system of nonlinear equations that express the polynomial exactness in order to 
find the nodes and weights; for instance, refer to [36]. Such methods have been found to yield 
good results for rather low-dimensional problems that feature symmetries and other invariance 
properties that can be exploited for simplification. The second class includes methods that 
search for polynomials that vanish at the nodes of the quadrature rule; for instance, refer 
to [20]. These methods have facilitated the study of important mathematical properties of 
quadrature rules using ideal and other theories. The third class consists of methods that involve 
the construction of quadrature rules by full or sparse-grid tensorization of suitable univariate 
rules; for instance, refer to [35]. These methods have found many successful applications in 
the context of stochastic expansion methods, but are limited in application to probability 
distributions that do not exhibit statistical dependence. Finally, methods belonging to the 
fourth class involve the selection of a subset of nodes from a given quadrature rule in order 
to construct a more efficient embedded rule. Many aspects of polynomial-based methods, 
including Theorems 5.1 and 5.2, are studied in [20, 29, 32, 34] and the references therein. 

Finally, it should be noted that references [37-40] propose alternative methods for the 
construction of embedded quadrature rules. The method given in [37, 38] relies on the 
computation of a basic optimal solution [26] to the linear program min ro ~^2 k=1 ro/c, subject 
to A-uj = b and vj > 0, where A and b are still defined by (38) and (39), respectively. A 
drawback of this method is that the existence of a basic optimal solution is generally not 
guaranteed. The existence of this solution is ensured [37] when the original quadrature rule 
has only positive weights and has the targeted degree of exactness, but a basic optimal solution 
may fail to exist in other cases because the imposed constraints may then be overly restrictive. 
The method given in [39, 40] relies on the solution of Avj = b for a basic solution using a 
QR factorization with column pivoting. This method is computationally less costly than the 
proposed method, but it can generally be expected to provide a quadrature rule for which the 
sum of the absolute values of the weights is larger and for which the error bound (40) is thus 
less favorable than the one for the quadrature rule obtained using the proposed method. 



6. Implementation 

In this section, we provide details on the implementation of the proposed dimension-reduction 
and measure-transformation methodology using stochastic expansion methods. 

6.1. Karhunen-Loeve decomposition 

Here, we describe the implementation of the KL decomposition of a random variable 
represented by a PC expansion and we show how this implementation in turn naturally provides 
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a PC expansion of the reduced random variables. Adopting the notations used in Sees. 2 and 3, 
we consider the construction of a reduced-dimensional representation of a second-order random 
variable q p with values in R w that is represented by a PC expansion of the following form: 

Q P = E 9«^«(€,C), 9«eK», (42) 

where {ip a , ot e N m+n } denotes a PC basis for the Hilbcrt space of Pi£ ^)-square-integrable 
functions from R m +™ into R with ipo — 1, as mentioned previously. We specifically consider the 
construction of the KL decomposition described in Sec. 3 that involves a weighting matrix W; 
the standard KL decomposition can be recovered easily by setting the weighting matrix equal 
to the identity matrix. Because of the orthonormality of the ip a , the mean vector q and the 
covariance matrix C q of q p can be deduced immediately from the PC coordinates as follows: 

q = Qo, (43) 
v 

C q = J2 *<*qI- (44) 

|c|=l 

Then, the solution of the generalized eigenproblem W T C q W& = \jW& provides the 
eigenvalues Ai > A 2 > . . . X w > and the associated eigenmodes 1 , . . . , (f> w required to 
construct a reduced-dimensional representation q p ' d of q p as follows: 

d 

q p4 =-q + J2^V^, (45) 
where the r] p are random variables with values in R such that 

v p = ^-(q p --qfw^ (46) 



and are zero-mean and uncorrelated. By substituting (42) in (46), a representation of each 
reduced random variable as a PC expansion is immediately obtained: 

p i 
V* = Vj,MtX) with Vj . a = —qlwtf, (47) 
|c|=i V X i 

thus indicating that the KL decomposition of a PC expansion naturally provides a complete 
probabilistic characterization of the reduced random variables as a PC expansion. 

It should be noted that the proposed methodology provides a representation of the solution 
and coupling variables associated with the subproblems as PC expansions in a combination of 
input random variables and reduced random variables of reduced-dimensional representations 
of exchanged information — refer to (8). Thus, when these reduced random variables are 
represented by PC expansions in the input random variables themselves — refer to (47), the 
proposed methodology requires the construction of KL decompositions of random variables 
that are represented by compositions of PC expansions. The implementation of such KL 
decompositions falls within the scope of the implementation mentioned above because the 
composition of two PC expansions can always be written equivalently as a PC expansion of 
the form (42), although as one that is truncated at a higher total degree. 
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6.2. Polynomial chaos with respect to arbitrary probability distributions 

Here, we provide details on the implementation of the Gram-Schmidt method described in 
Sec. 4. Adopting the notations used in Sec. 4, we consider the construction of a set of z-variate 
orthonormal polynomials {1^,0 < |/3| < q} up to a specific total degree q with respect to the 
probability distribution P x on the finite-dimensional Euclidean space R 2 . To ensure that the 
Gram-Schmidt method can be executed properly, we assume that P x satisfies the conditions 
given by (25) and (26). Let the z-variate monomials {x 13 , < |/3| < q} then be ordered in a 
sequence x^ 1 > • ■ • > X^ 1 * i where fj, = dim(7 , |) denotes the number of monomials in this sequence, 
as mentioned previously. Following the approach given in [22, 41], we propose to implement the 
Gram-Schmidt orthonormalization of this sequence using a Cholesky factorization of its Gram 
matrix. The Gram matrix G of the sequence of monomials x^ 1 > • • • > is the ^.-dimensional, 
square, and symmetric matrix that collects the inner products of these monomials as follows: 



G 



J R ,X^X^dP x ... / R , X ft X^dP x 
J^X^X^dP x ... t^^^. 



(48) 



Owing to the assumption that the conditions given by (25) and (26) are fulfilled, the Gram 
matrix G is bounded and positive definite. The Cholesky factorization G — R T R of the Gram 
matrix thus provides a ^-dimensional, square, and upper triangular matrix R with strictly 
positive diagonal entries. The inversion of this upper triangular matrix R in turn provides 
an upper triangular matrix RT 1 . Let Sij denote the entries of then, the orthonormal 

polynomials {Tp,0 < \/3\ < q} to be determined are obtained as follows: 

(x) = si 3 X 01 + S2jX^ 2 + ■■■ + SjjX 0i , 3 = 1, • • • , W (49) 

in fact, with Xu = [x^ 1 ■ ■ ■ X M ] T an d T M = [Tp 1 . . . r i g^] T , we immediately obtain 

/ T^ X )T^x) T dP x - R T f X^xldP x R 1 = R T GR ' = I. (50) 

JR* JR* 

Here, I denotes the ^.-dimensional identity matrix. 

It should be noted that the computation of the basis of multivariate polynomials 
by orthonormalization of the basis of monomials may be overly sensitive to numerical 
approximation and roundoff errors for large values of the dimension z and total degree q. 
This issue can be addressed, for instance, by following the approaches given in [22, 42]. 

6.3. Quadrature rules with respect to arbitrary probability distributions 

Here, we provide details on the method for the computational construction of embedded 
quadrature rules described in Sec. 5. At the core of this method is the solution of the L 1 - 
minimization problem (37) for a sparse optimal solution. Optimization theory provides two 
types of approaches for solving L 1 -minimization problems. Both types of approaches first 
involve the reformulation of the i 1 -minimization problem as an equivalent linear-programming 
problem; then, the approaches of the first type apply a simplex algorithm, whereas the 
approaches of the second type apply an interior-point algorithm [26, 27]. The application 
of a simplex algorithm directly provides a sparse optimal solution because a simplex algorithm 
computes an optimal solution by performing a sequence of pivoting operations on basic feasible 
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solutions. In contrast, the application of an interior-point algorithm generally provides a fully 
populated optimal solution; thus, the application of an interior-point algorithm requires an 
additional computational procedure to extract a sparse optimal solution [27, 28]. In this section, 
we provide details on the solution of (37) using an interior-point algorithm. It should be noted 
that our discussion can readily be extended to obtain a solution to (37) by a simplex algorithm; 
however, for brevity, we do not explicitly demonstrate this extension. 

Adopting the notations used in Sec. 5, we consider the construction of a quadrature rule for 
integration with respect to the probability distribution P x on the finite-dimensional Euclidean 
space R z with accuracy level A. To ensure that the proposed method can be executed properly, 
we assume that the probability distribution P x satisfies the condition given by (33). 

Let us assume that a quadrature rule for integration with respect to P x is already available 
and this quadrature rule has a very large number (denoted by v) of nodes and associated 
weights. Then, we focus on the construction of an embedded quadrature rule with fewer 
nodes and weights using the solution of the ^-minimization problem min OTeR s ||^|| £ i subject 
to Avj = 6, in which A is the (i x i>-dimensional matrix and 6 the /i-dimensional vector 
defined by (38) and (39), respectively, introduced to impose the prescribed accuracy level A; it 
should be noted that \i = dim('P:? A-1 ) denotes the number of imposed equality constraints, as 
mentioned previously. Following the approach given in [26, 27], we convert the ^-minimization 
problem into a linear program as follows: 

min e T i, subject to Avj = b, vj — t < 0, and vj + t > 0, (51) 

where t is a ^-dimensional vector of slack variables and e is a ^-dimensional vector defined 
by e = [1,...,1] T . The Karush-Kuhn- Tucker optimality conditions for the Lagrangian 
associated with the constrained optimization problem (51) are then expressed as follows: 

A T \ — n + fi = and e — /j, — fi = (stationarity), (52) 

Avj 6 = 0, —z& + t>0, and vj + 1 > (primal feasibility), (53) 

(i > and fi > (dual feasibility), (54) 

^ife(— tufe + tk) =0 and fiki^k + ifc) = 0, 1 < k < D (complementarity), (55) 

where the /i-dimensional vector A and the ^-dimensional vectors /x and fi collect the Lagrange 
multipliers. After eliminating fi = e — /x, these optimality conditions are expressed as follows: 

A T \ 2n + e = (stationarity), (56) 

Azu -6 = 0, + 1 > 0, and za + 1 > (primal feasibility), (57) 

fi > and e — fi > (dual feasibility) , (58) 

Hk(— tJk + tk) = and (1 — Hk){^k + tk) = 0, 1 < k < v (complementarity). (59) 

The primal linear program (51) is associated with the dual linear program 

min 6 T A, subject to A T \ — 2/x + e = 0, /x > 0, and e /x > 0; (60) 

the relationship between the primal linear program (51) and the dual linear program (60) 
is that their Karush-Kuhn- Tucker optimality conditions are identical. Consequently, using 
the complementarity (59), which implies that /x T (— w + t) + (e — fi) T (zv + t) =0, the 
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stationarity (56), which implies that —X t Avj + e T t = 0, and the primal feasibility (57), 
it can be shown that the optimal values for the primal and dual linear programs are the same: 

A T b = e T i. (61) 

Further, it should be noted that the complementarity (59) implies that fx k — when w k < 
and that fj, k = 1 when vj k > 0; the values of the Lagrange multipliers [i k associated with the 
nonzero components m k are thus determined solely by the sign of these components. 

Now, the solution of the linear program using a primal-dual interior-point algorithm yields 
a primal-dual optimal solution (zu,t,X,/j,) that solves the Karush-Kuhn- Tucker optimality 
conditions (56)-(59). However, the vector vj obtained as a part of such a quadruple, 
i.e. (vj, t, A, fi), is in general fully populated and therefore not sparse. 

Following the approach given in [26, 27], we propose to use the procedure described next to 
extract a sparse optimal solution. Let £> U N be a partitioning of the index set {1, . . . , v} such 
that B collects the indices associated with the nonzero components of vj and M collects the 
remaining indices associated with the vanishing components of vj. Let the number of nonzero 
components be greater than the number (denoted by fi) of imposed equality constraints. Then, 
there exists necessarily a vector z ^ such that Az = and the components of z labeled by 
the indices in N vanish; and there exists necessarily a scalar a such that vj + az has more 
vanishing components than vj and each component w k + az k either vanishes or has the same 
sign as w k . Because the components m k + az k either vanish or have the same sign as the 
components vjk, the Lagrange multipliers fi k satisfy 

ti k (-(-cu k +az k ) + \m k +az k \j = and (1- ^ k )(^(zu k +az k ) + \w k +az k \j =0, 1 < k < D; (62) 

hence, the quadruple (w + az, \w + az\, A, /x) is also a primal-dual optimal solution that 
satisfies (56)-(59). Consequently, the quadruple (w + az, \w + az\, A, /i) satisfies — X t A(vj + 
az) + e T \vj + az\ = such that, with Az — 0, the value taken by the objective function 
at (vj + az, \vj + az\, A, fi) is equal to the value taken by the objective function at (tu, t, A, fx): 

e T \zu + az\ = X t A(vj + az) = A T Azv = X T b = e T t = e T \zv\. (63) 

In conclusion, this procedure yields a vector vj + az that at least has one nonzero component 
less than vj but is still optimal in that the value of the objective function remains unchanged. 
Clearly, the repeated application of this procedure will yield an optimal solution that has at 
most as many nonzero components as there are imposed equality constraints. Algorithm 1 
outlines the method thus obtained for the computation of an embedded quadrature rule. 



7. Realization for a multiphysics problem 

7.1. Problem formulation 

[Figure 1 about here.] 

We consider the stationary transport of neutrons in a one-dimensional reactor with 
temperature feedback [43]. Let the reactor occupy an open interval ]0, L[ (Fig. 1). The problem 
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Input : Quadrature level A; 

Quadrature rule {(x k ,Wk), 1 < fc < i>} for integration with respect to P x 

^-minimization 

Select polynomial basis 7Ti, . . . , 7r M for "Pf A_1 with fj, = dim^^ -1 ); 
Construct matrix A using (38) as follows: 



A = 



.^(xi) ^(x 2 ) • 

Construct vector b using (39) as follows: 

b=[J Rz n 1 (x)dP x . 
Use interior-point algorithm to solve linear program 

subject to [A T -21] 



Tl(Xi>-l) TTlfe)' 

^(Xp-i) ^(x&). 



/ R . *n(x)dP x ] T ; 



min [-6 Ol 

AeR^GR 5 



= — e and < fi < e, 



for primal-dual optimal solution (vo, t, A, ft); 



end 

Extraction of sparse optimal solution 
repeat 

Partition {1, . . . , i>} as B U N such that w k = if k is in A/"; 
Find vector z such that Az = and z k — if k is in A/"; 
Set a = min{|z fe /n7A;| : k e B , z k ^ 0, and sign(z fe ) 7^ sign(tu fe )}; 
Update vj to + az; 
until a = 0; 
end 

Construction of embedded quadrature rule 

Synthesize quadrature rule {{xt = Xk,i w l = ro fcf)i 1 — ^ ^ 
in which £> = {fci, . . . , fc„} with necessarily v < [i = dim(7 , ^ A_1 ); 
end 

Algorithm 1: Computation of an embedded quadrature rule. 



then involves finding the temperature T and neutron flux $ such that 

I K>f ) - 

under homogeneous Neumann boundary conditions. The first term on the left-hand side of the 
heat subproblem represents heat conduction, and the second term represents the transmission 
of heat to the surroundings; further, the right-hand side represents a distributed heat source 



-£ f £ f (T)$, 
i/E f (T))$ = -a, 



(64) 
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proportional to the neutron flux. The first term on the left-hand side of the neutronics 
subproblem represents neutron diffusion, and the second term represents the net effect of the 
absorption and generation of neutrons; further, the right-hand side represents a distributed 
neutron source. The coefficients k and h are the heat conductivity and heat transmittivity, 
respectively; the temperature T m is the ambient temperature; and v and Ef are the number 
of neutrons and the energy released per fission reaction, respectively. The coefficients D, S a , 
and Sf are the neutron diffusion constant, fission cross section, and absorption cross section, 
respectively; these coefficients depend on the reactor temperature as follows: 



D(T(x)) = D rel] J^, S a (T(x)) = S a , rcf y|^, E f (T(z)) = (65) 
7.2. Deterministic weak formulation 

Let H = H 1 (]0,L[) be the space of functions that are sufficiently regular to describe the 
solutions to the heat and neutronics subproblems. The weak formulation then involves finding T 
and $ in H such that 



/ k——dx+ h{T-T ca )Sdx = E t Y, t {T)^Sdx, V5 € H 1 

Jo dx ax J Q J 

/ D(T)- — —dx+ [E a (T)-i/Ef(T))$*da;= / s^dx, V* e H. 

Jo dx dx Jo V ) J 



(66) 



7.3. Random thermal transmittivity 

We incorporate uncertainties by modeling the thermal transmittivity as a random 
field {h(x, ■), 1 < x < L} such that 



(m 



(67) 



where the £j are statistically independent uniform random variables defined on a probability 
triple (9, T, P) with values in [— f , f ] and the V3^j are thus uniform random variables with unit 
standard deviation; further, the Xj and <j>> are the eigenvalues and eigenmodes, respectively, 
of the eigenproblem C^) — Xjfl , where C is the covariance integral operator with 

^/ \ 4a 2 . 9 f ir(x — y)\ 

C ^ = ^W Sm ( A 2^ Z J (68) 

as the kernel; here, the parameter a is the spatial correlation length of {h(x, ■), 1 < x < L}. 
Clearly, the random field {h(x, -),1 < x < L} thus obtained is such that the random 
variable h(x, •) has the mean h and coefficient of variation S at every position x, at least 
when the approximation error introduced owing to the truncation of the expansion after m 
terms is not taken into account. 
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7.4- Stochastic weak formulation 

The weak formulation of the stochastic problem involves finding random variables T and <J> 
defined on (0,T, P), with values in H, such that 



/ k ^TlT dx + / K£){ T - T o°)Sdx= / EfEf(T)$Sdx, VS G H, 
Jo dx dx J Q J Q 

r L r/$ d^ f L / \ r L 

/ D(T)- — —dx+ [E a (T)-i/Ef(T))$*da;= / s*dx, V* e £T. 
Jo « x dx J Q \ J J 



(69) 



7. 5. Discretization of space 

The finite element (FE) method is used for the discretization of space. The domain [0, L] is 
meshed using r — 1 elements of equal length. Let N\,. . . ,N r then be a basis of element- wise 
linear shape functions such that Nj takes value 1 at the j-th node and at other nodes. Using 
this basis, the random temperature T and neutron flux $ are approximated as follows: 

r 

7 1 (70) 

The FE discretization of the stochastic weak formulation (69) then involves finding random 
vectors T = {Ti, . . . , T r } and $ = {$i, . . . , $ r } defined on (0, T, P), with values in M r , which 
collect the nodal values of the random temperature and neutron flux such that 

[K + H{$]T = q(*,T), 

[D(T) + M(T)]* = s. [ ' 

Here, K, H, D(T), and M(T) are r-dimensional matrices, and </(<&, T) and s are r- 
dimcnsional vectors such that 

SjKS 2 = [ L k^^dx, (72) 
J Q dx dx 

SjHS 2 = f hS[S r 2 dx, (73) 
Jo 

•fD<m.-jfi>Cr)££M.. (w) 

*JM(T)* 2 = / (s.fT'-)-^,^'-))*;*^, (75) 



S T q(T,$>) = EiY,{(T r )<$> r S r dx+ hT^STda 
Jo Jo 



(76) 



S T s= / «tf r da;. (77) 
Jo 
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7.6. Reformulation as a realization of the model problem 

The aforementioned illustration problem can be reformulated as a particular realization of the 
general model problem introduced in Sec. 2 as follows: 



r = a(T,*,£), a: 
* = b{T), b : ] 



(78) 



where a(T,*,£) = [X>ff(£)] _1 g(*,T) and b(T) = [D(T)+M(T)]~ 1 s. This reformulation 
indicates that the illustration problem is a simplified realization of the model problem, for three 
reasons. First, the data of the neutronics subproblem are not affected by their own sources 
of uncertainty Second, the neutronics subproblem admits a direct solution that does not 
require iteration. Lastly, the neutronics and heat subproblcms are coupled directly through 
their solution variables rather than through intermediate coupling variables. 

7.7. Dimension reduction 

Now, we will demonstrate the proposed methodology by approximating the random 
temperature by a truncated KL decomposition as it is communicated from the heat to the 
neutronics subproblem. At iteration I, let the random temperature be represented by a PC 
expansion as follows: 



The mean and covariance of T l ' v are then given by 



re 



P 



|a|=l 



(79) 

(80) 
(81) 



Further, let the r-dimcnsional square matrix W be the Gram matrix of the FE basis, i.e., 



w 



(N r , N- 



1 H 



(N r ,N, 



r H 



(82) 



where the inner product (•, is such that (Si, S2)h = Jq SiS^cfcc + f$ (dSi / dx)(dS2 / dx)dx 
for any pair Si and S2 of functions in H. The solution of the generalized eigenproblem 

C^pW <j) : '' l = XjWcfy''^ then provides the eigenvalues and the associated eigenmodes 
required to construct a reduced-dimensional representation T £ 'P' d of T e ' p as follows: 



where the r{- v are random variables defined on (9, T, P), with values in M, such that 



Vj" 



1 =(T'^r)Vf { . 



(83) 



(84) 
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By substituting (79) in (84), a representation of the r\\' v as a PC expansion is obtained: 



J* = £ < 

H=i 



Vv*(£) with rf c 



(5*c 



(85) 



thus completely characterizing the reduced random variables as a PC expansion. 

It should be noted that the random neutron flux, in principle, could also be reduced as it 
passes from the neutronics subproblem to the heat subproblem. However, because the data of 
the neutronics subproblem are not affected by their own sources of uncertainty, a reduction of 
the random neutron flux would not lower the number of sources of uncertainty that enter the 
heat subproblem and thus would not lead to a solution of the heat subproblem in a reduced- 
dimensional space. This extension is therefore not demonstrated. 



7.i?. Measure transformation 

Whereas the random variables £ = (£i , . . . , £ m ) necessarily constitute the sources of uncertainty 
that enter the heat subproblem, the reduced random variables rf — {r)[ , . . . ,7]^) of the KL 
decomposition of the random temperature can be construed as the sources of uncertainty that 
enter the neutronics subproblem. Then, the proposed methodology leads to the approximation 
of the random temperature and neutron flux by PC expansions as follows: 

\a\=0 
|/3|=0 

i.e., we obtain the approximation of the random temperature by a PC expansion in the input 
random variables and the approximation of the random neutron flux by a PC expansion in 
the reduced random variables. 

We select the polynomials {V'a, < |cc| < p} as normalized Legendre polynomials, and 
we construct the polynomials {r^,0 < |/3| < q} at each iteration using the method given 
in Sees. 4 and 6.2. Further, we select the quadrature rule for integration with respect to the 
probability distribution of the input random variables, denoted as {(£ k , v k ), 1 < k < N}, as a 
sparse-grid Gauss-Legendre quadrature rule of level p+1, and we construct the quadrature rule 
for integration with respect to the probability distribution of the reduced random variables, 
denoted as {(r] l k , w k ), 1 < k < j/}, using the method given in Sees. 5 and 6.3 for level A = q+1. 

It should be noted that because the polynomial transformation (85) is continuous and 
the probability distribution of the input random variables has a closed and bounded 
support [—1,1]™, the probability distribution of the reduced random variables also has a closed 
and bounded support. With reference to Sees. 4.4 and 5.5, this property suffices to ensure the 
convergence of the PC expansion of the random neutron flux as the total degree q increases. 



(86) 



7.9. Implementation 

The abovementioned computational construction of the polynomials and quadrature rules 
requires that integrals be computed with respect to the probability distribution of the 
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reduced random variables to obtain the requisite Gram matrix and moment vector. Because 
the probability distribution of the reduced random variables is characterized by the KL 
decomposition (83) as a transformation of the probability distribution of the input random 
variables through the PC expansion (85), these integrals can be evaluated easily using a 
quadrature rule for integration with respect to the probability distribution of the input random 
variables after a "change of variables." 

Owing to the PC expansion (85), the integral of a function / from R d into R with respect 
to the probability distribution P^ of the reduced random variables can be reformulated as an 
integral with respect to the probability distribution P$ by a "change of variables" as follows: 

/ f(v)dPt= I fU^tjjdPs, (87) 

provided that either integral is well defined [44]; thus, the numerical evaluation of this integral 
can be performed using a quadrature rule for integration with respect to 

In this work, we use the Monte Carlo method to numerically evaluate the entries of the 
Gram matrix and the components of the moment vector as follows: 

- , MC 

„ 1 MC 

b e = ^ tf*f idP „ _J_ J2 (^ k )f\ (89) 

where {£ fe , 1 < k < MC} collects MC independent samples of the input random variables. 

In addition, the method given in Sees. 5 and 6.3 requires a quadrature rule to provide a 
collection of candidate nodes from which the nodes of the embedded quadrature rule can be 
selected. We adopt the following construction in this work. First, we construct a sparse-grid 
Gauss-Legendre quadrature rule, denoted as {(£&, Wk), 1 < k < z>}, for integration with respect 
to the probability distribution of the input random variables. Next, we carry out a "change 
of variables" to transform this sparse-grid Gauss-Legendre quadrature rule to a corresponding 
quadrature rule of the form {(t)|, u>k), 1 < k < v} for integration with respect to the probability 
distribution of the reduced random variables by choosing the nodes as rf k — rj e ' p (£ k ) and 
keeping the weights unchanged. Finally, from this quadrature rule, we construct an embedded 
quadrature rule {(ri^w^), 1 < k < v} of the desired level A by applying Algorithm 1. 

The aforementioned sparse-grid Gauss-Legendre quadrature rule is parameterized itself by 
its own level, denoted here as A. As a larger A is chosen, the number of nodes that the sparse- 
grid Gauss-Legendre quadrature rule has increases; thus, an increase in A provides Algorithm 1 
with a greater choice of candidate nodes to select from, and therefore, Algorithm 1 can be 
expected to yield a better embedded quadrature rule that provides a smaller value for the sum 
of the absolute values of the weights. Because this sum has a lower bound of 1, the adequate 
level A can readily be selected by monitoring the convergence of this sum with respect to A. 

This is not the only implementation available. The Monte Carlo method can be replaced 
by a fully tensorized or sparse-grid quadrature rule to evaluate (88) and (89). Further, the 
Gram matrix and moment vector can be deduced using polynomial algebra from the PC 
expansion (85) and the moments of P$. Finally, the sparse-grid Gauss-Legendre quadrature 
rule can be replaced by Clenshaw-Curtis, Monte Carlo, or other rules. 
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7.10. Selection of the reduced dimension and the polynomial degree 

At each iteration, we select the number of terms retained in (83) by the KL 
decomposition T l - p - d of T l - p as the smallest dimension d that satisfies the following condition: 



w 



dP < en 



dP, 



We 



w 



(90) 



where ei is a prescribed tolerance level. Further, at each iteration, we truncate the PC 
expansion <I? £ ' P in (86) at the smallest total degree q that satisfies the following condition: 



&t,q _ $l,g-l 



W 



dP < e 2 \ 



w 



dP, 



W e N, 



(91) 



where £2 is a prescribed tolerance level. Clearly, these criteria may result in the dependence of 
the reduced dimension d and the total degree q on the iteration i. 



7.11. Concluding remarks 

Algorithm 2 summarizes the implementation of the illustration problem presented in this 
section. The key feature of this implementation is that it enables a solution of the neutronics 
subproblem in a reduced-dimensional space when the KL decomposition can extract a low- 
dimensional representation of the random temperature (d < m), while maintaining accuracy. 

The solution of the neutronics subproblem in a reduced-dimensional space can be expected 
to reduce the number of terms that are required in the PC expansion of the random neutron 
flux to achieve sufficient accuracy. Further, the solution of the neutronics subproblem in a 
reduced-dimensional space can be expected to reduce the number of quadrature nodes that 
are required for the nonintrusive projection method to achieve sufficient accuracy in the PC 
coordinates of the random neutron flux and therefore to reduce the number of times a sample 
of the neutronics subproblem must be solved, thus lowering the computational cost. 



8. Numerical results 

We obtained numerical results using the following properties. We assumed the reactor 
to have a length of L = 100 [cm]. Further, we assumed a deterministic and position- 
independent heat conductivity k — 100 [J/K/cm/s]; ambient temperature = 390 [K]; fission 
energy Ef = 3.0.E-11 [J/neutrons]; fission cross section £ a ,rcf = 0.0075 [cm -1 ]; neutron-diffusion 
constant D re f = 2.2 [cm]; absorption cross section £ a ,ref = 0.0195 [cm -1 ]; multiplication 
factor v = 2.2; neutron source s = 5.0iTl [neutrons/s/cm 3 ]; and temperatures T ro f = 390 [K], 
T m in = 390 [K], and T max = 1000 [K]. 

[Figure 2 about here.] 

In addition, we used a thermal transmittivity random field with position-independent 
mean h — 0.17 [J/K/cm 3 /s], spatial correlation length a — 15 [cm], and coefficient of 
variation 5 = 10%. We retained to = 10 terms in expansion (67). Figure 2(a) shows a few 
sample paths of the random field {h(x, •), < x < L} thus obtained. Figure 2(b) shows the 
10 largest magnitude eigenvalues of the covariance integral operator. 
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Input : Error tolerance levels ei and 62; 

PC basis {ipa, < \a\ < p) up to total degree p w.r.t. Pj; 
Quadrature rule {(£ k ,v k ), 1 < k < N} of level p + 1 w.r.t. P^; 

repeat 

heat subproblem 

for k = 1 to A^+i do 

Solve [K + Jf(* fc )]T%) =g(T^ 1 ^(^),^- 1 ^(4)), 

with s'- 1 "^) = Er/aNo^-^-^EfaNi^aVa^)); 
end 

Compute PC coordinates of using T' a = E^I* r'^V'e. 
end 

dimension reduction 

Compute mean T 1 = T l and covariance matrix C~ = Eft* 1=1 T a(T e a ) T ; 
Solve eigenproblem W T C^Wft* 1 = X^Wcj) 3 ^; 
Choose d such that ^/E 3 C +1 < £i\/Ef a | =1 (^) T W^; 
Compute PC coordinates of r[- v by = (T^,) T W0 J "'* for j = 1 to rf; 
end 

neutronics subproblem 

9 = 0; 
repeat 

measure transformation 

Compute PC basis {Pa, < |/3| < 17} up to total degree q w.r.t. P^; 
Compute quadrature rule {(77^, 1 < fc < ;/} of level q + 1 w.r.t. P^; 
end 

for fc = 1 to 1/ do 

Solve [D(T^ d (r,i)) + M(T*-*V*))]*V*) = a , 
with T*-» V*) = r £ + E U \f^rf, k <P' l \ 
end 

Compute PC coordinates of $ e ' q using ^ = Efeli * < ( T 7fc)rjs('7l) u; k> 
9 = 9+1; 

until (^E|^| =q II^Hw<^\/E^| =0 ll^ll 2 w); 
end 

^ = *+l; 
until (convergence); 

Algorithm 2: Implementation of the illustration problem. 
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8.1. Monte Carlo sampling implementation 

[Figure 3 about here.] 

First, we carried out a Monte Carlo simulation. We generated MC — 100, 000 sample 
paths of the thermal transmittivity random field. Then, for each of these sample paths, we 
constructed the associated deterministic multiphysics model, each of which we solved using the 
FE method for the spatial discretization and Gauss-Seidel iteration as the iterative method. 
We systematically obtained converged results for r — 1 = 40 finite elements and 20 iterations. 

Figure 3 shows a few samples of the random temperature and neutron flux thus obtained. 
We can observe that the samples of the random temperature (Fig. 3(b)) are smoother than the 
samples of the thermal transmittivity random field (Fig. 2(a)); i.e., the former exhibit less rapid 
oscillations with respect to the position in the reactor than the latter. We had demonstrated 
in [11] that this behavior can be attributed to the large magnitude of the diffusion term of the 
heat subproblem which reduces the nonuniformity of the samples of the random temperature. 

8.2. PC-based implementation involving dimension reduction and measure transformation 

Next, we implemented the proposed PC-based iterative method involving dimension reduction 
and measure transformation. This implementation exactly corresponded to Algorithm 2. We 
obtained the results to follow by systematically setting the total degree of the PC expansion 
of the random temperature to p = 4 and, with reference to (90) and (91), using a range of 
values for the error tolerance levels e\ and e 2 adopted to determine the reduced dimension and 
the total degree of the PC expansion of the random neutron flux at each iteration. We discuss 
convergence as a function of these error tolerance levels later. Now, we present detailed results 
obtained for ei = 0.05 and e 2 = 0.000001. 

[Figure 4 about here.] 

Figure 4 shows the convergence of the iterative method as a function of the number of 
iterations; note that the superscript oo is used in the figure captions to indicate convergence 
with respect to the number of iterations. The iterative method converged at a linear rate up 
to approximately iteration I = 10, after which linear-solver tolerances became dominant and 
prevented further convergence. All results to follow were obtained at iteration I = 20 and can 
thus be considered to have converged with respect to the number of iterations. 

[Figure 5 about here.] 

[Figure 6 about here.] 

Figure 5 shows a few components of the KL decomposition of the random temperature, 
namely, the mean temperature, the 10 largest magnitude eigenvalues, and the eigenmodes 
and reduced random variables associated with the two largest magnitude eigenvalues. We can 
observe that the eigenvalues of the KL decomposition of the random temperature (Fig. 5(b)) 
decay at a faster rate than those of the KL decomposition of the thermal transmittivity random 
field (Fig. 2(b)), which is consistent with our earlier observation that the samples of the former 
are smoother than those of the latter. 

Figure 6 shows the joint and marginal probability density functions of the reduced random 
variables. Clearly, the joint probability density function exhibits statistical dependence and 
the marginal probability density functions are not "labeled." 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



28 



M. ARNST, R. GHANEM, E. PHIPPS AND J. RED-HORSE 



[Figure 7 about here.] 
[Figure 8 about here.] 

At iteration £ = 20, a KL decomposition retaining only d = 2 terms was found to be 
sufficiently accurate to satisfy (90) for ei = 0.05; thus, at iteration I = 20, the measure 
transformation necessitated the construction of orthonormal polynomials and quadrature rules 
with respect to the probability distribution of the first and second reduced random variables. 

Figures 7 and 8 illustrate the proposed computational construction of orthonormal 
polynomials. With reference to Sec. 7.9, Fig. 7 demonstrates the convergence of the Monte 
Carlo estimate of an entry of the requisite Gram matrix with respect to the number of samples. 
Figure 8 shows the obtained orthonormal polynomials up to a total degree of q = 3. 

[Figure 9 about here.] 

[Figure 10 about here.] 

[Figure 11 about here.] 

[Figure 12 about here.] 

Figures 9-12 illustrate the proposed computational construction of quadrature rules. With 
reference to Sec. 7.9, Fig. 9 demonstrates the convergence of the Monte Carlo estimate of a 
component of the requisite moment vector with respect to the number of samples. 

We obtained the results shown in Fig. 10 as follows. First, we constructed three sparse- 
grid Gauss-Legendre quadrature rules of the form {(£ fc , Wk), 1 < k < v} for integration with 
respect to the probability distribution of the input random variables: the first rule was of the 
level A = 3 and had v = 261 nodes; the second was of the level A = 4 and had v — 2, 441 nodes; 
and the third was of the level A = 5 and had v = 18, 881 nodes. Next, we carried out a "change 
of variables" to transform each of these sparse-grid Gauss-Legendre quadrature rules into a 
corresponding quadrature rule of the form {(rf k — »/' p (£ fc ), il^), 1 < k < v} for integration 
with respect to the probability distribution of the reduced random variables; Fig. 10((a), (d), 
and(g)) shows the nodes of the quadrature rules thus obtained. Finally, on the basis of each of 
these quadrature rules, we constructed an embedded quadrature rule {(ij e k7 Wk), 1 < k < v 1 } 
of level A = 4 using Algorithm 1; Fig. 10((b), (e), and(h)) and 10((c), (f), and(i)) shows the 
nodes and weights of the embedded quadrature rules thus obtained. Figure 10((c), (f), and(i)) 
indicates that as the level A was increased, Algorithm 1 was provided with a greater choice 
of candidate nodes from which the nodes of the embedded quadrature rule could be selected 
and thus provided a better embedded quadrature rule with a smaller value for the sum of 
the absolute values of the weights. Figure 11 demonstrates the convergence of the sum of the 
absolute values of the weights with respect to the level A, indicating convergence at A = 5. 

Figure 12 shows the nodes and weights of the embedded quadrature rules of levels A = 1, 
2, and 3 constructed similarly on the basis of the sparse-grid Gauss-Legendre quadrature rule 
of level A = 5. This figure indicates that as the level A was increased, higher accuracy was 
required, and thus an embedded quadrature rule with more nodes was systematically obtained. 

[Figure 13 about here.] 
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At iteration I — 20, a PC expansion truncated at a total degree of q = 3 was sufficiently 
accurate to satisfy (91) for e 2 = 0.000001. The representation of the random temperature by a 
PC expansion of dimension m = 10 and total degree p = 4 requires 1,001 = 14!/10!/4! terms, 
whereas the representation of the random neutron flux by a PC expansion of dimension d = 2 
and total degree q = 3 requires only 10 = 5!/2!/3! terms. Figure 13 shows a few PC coordinates. 

[Figure 14 about here.] 

Figure 14 shows a few samples of the random temperature and neutron flux deduced from 
the PC expansions obtained as the output of the solution algorithm. The samples of the input 
random variables used to synthesize the samples of the random temperature and neutron flux 
shown in Fig. 14 were identical to those used to generate the samples shown in Fig. 3. The 
similarity of the samples in Figs. 3 and 14 indicates that the PC-based surrogate model not 
only provides an accurate global representation of the multiphysics model but also is capable 
of accurately reproducing a sample-wise response. 

8.3. Convergence analysis 

[Figure 15 about here.] 
[Figure 16 about here.] 

We repeated the PC-based simulation for several values of the error tolerance levels; each 
of these values corresponded to specific accuracies that the KL decomposition of the random 
temperature and the PC expansion of the random neutron flux were required to maintain 
at each iteration. Figures 15(a) and 16(a) indicate that this KL decomposition and the PC 
expansion systematically retained more terms when higher accuracy was required. 

Further, Figs. 15((b) and (c)) and 16((b) and (c)) indicate that the distance between the 
successive approximations determined by the Monte Carlo and PC-based iterative methods 
remained bounded as the iterations progressed and this distance can be reduced systematically 
by improving the accuracy of the KL decomposition of the random temperature and the PC 
expansion of the random neutron flux by decreasing the respective error tolerance levels. 

8.4- Concluding remarks 

The proposed methodology enabled the solution of the neutronics subproblem in a reduced- 
dimensional space because the KL decomposition was able to extract a low-dimensional 
representation of the random temperature as it passed from the heat subproblem to the 
neutronics subproblem. While accuracy was maintained, the solution of the neutronics 
subproblem in a reduced-dimensional space resulted in computational gains, for two reasons. 
First, it enabled the accurate representation of the random neutron flux by a PC expansion 
that had only a few terms. Second, it enabled the computation of the coordinates in the PC 
expansion of the random neutron flux using a quadrature rule that had only a few nodes, thus 
requiring the solution of only a few samples of the neutronics subproblem per iteration. 
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9. Conclusion 

While most coupled models can be expected to be affected by a large number of sources 
of uncertainty, information exchanged between subproblems and iterations often resides in a 
considerably lower dimensional space than the sources themselves. In this work, we thus used 
a dimension-reduction technique to extract a low-dimensional representation of information as 
it passes from subproblcm to subproblem and from iteration to iteration, and we proposed 
measure-transformation techniques that allows implementations to exploit this dimension 
reduction to achieve a computationally efficient solution of the subproblems in a reduced- 
dimensional space. The effectiveness of the proposed methodology was demonstrated on a 
multiphysics problem relevant to nuclear reactors. 
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Figure 1. Schematic representation of the problem. 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cls 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



36 



FIGURES 



0.22 




Position [cm] 

(a) Samples. 

40 1 1 , 



30- 




2 4 6 8 10 

Index [-] 

(b) Eigenvalues. 

Figure 2. Thermal transmittivity random field: (a) five samples and (b) ten largest magnitude 

eigenvalues of the covariance integral operator. 
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Figure 3. Monte Carlo simulation: five samples of the solution. 
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Figure 5. PC-based simulation: mean, eigenvalues, first and second eigenmode, and PC coordinates of 
the first and second reduced random variables of the KL decomposition of the random temperature. 
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Figure 6. PC-based simulation: joint and marginal probability density functions of the first and second 
reduced random variables of the KL decomposition of the random temperature. 
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Figure 8. PC-based simulation: computed orthonormal polynomials up to a total degree of q = 3. 
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Figure 10. PC-based simulation: nodes of the quadrature rules {(iji = f7 f,p (£fc), Wk), 1 < k < v} and 
nodes and weights of the embedded quadrature rules {(rj k ,Wk), 1 < k < u } of level A = 4 obtained 
from sparse-grid Gauss-Legendre quadrature rules {(£ fe , u>k), 1 < k < i>} of levels A = 3, 4, and 5. 
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Figure 12. PC-based simulation: embedded quadrature rules {{jf k , Wk), 1 < < of levels A = 1, 2, 
and 3 obtained from the sparse-grid Gauss-Legendre quadrature rule of level A = 5. 
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Figure 13. PC-based simulation: PC coordinates at x = 10 [cm] of the solution. 
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Figure 14. PC-based simulation: five samples of the solution. 
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Figure 15. Convergence analysis: (a) reduced dimension; and (b) and (c) mean-square distance between 
the successive approximations determined by the Monte Carlo and PC-based iterative methods for 
ei = 0.20 (circles), ei = 0.05 (squares), and ei = 0.01 (diamonds) and 62 = 0.0001 as a function of the 

iteration. 
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Figure 16. Convergence analysis: (a) total degree; and (b) and (c) mean-square distance between the 
successive approximations determined by the Monte Carlo and PC-based iterative methods for ei = 
0.05 and £2 = 1 (circles), £2 = 0.01 (squares), £2 = 0.0001 (diamonds), and £2 = 0.000001 (triangles) 

as a function of the iteration. 
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